home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / rexx / imc / rexx-imc.5 / rxmathfn.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-06-25  |  7.9 KB  |  299 lines
  1. /* A sample function package for REXX/imc to provide the math functions */
  2. /*                                                 (C) Ian Collier 1992 */
  3.  
  4. #include"functions.h"
  5. #include<math.h>
  6. #include<stdlib.h>
  7.  
  8. /* Below is the dictionary which informs REXX what functions are available
  9.    and where to find them */
  10.  
  11. int rxsin(),rxcos(),rxtan(),rxexp(),rxatn(),rxln(),rxpower(),rxasn(),
  12.     rxacs(),rxsqr();
  13.  
  14. dictionary rxdictionary[]=
  15. {  "ACOS"   ,rxacs,
  16.    "ASIN"   ,rxasn,
  17.    "ATAN"   ,rxatn,
  18.    "COS"    ,rxcos,
  19.    "EXP"    ,rxexp,
  20.    "LN"     ,rxln,
  21.    "SIN"    ,rxsin,
  22.    "SQRT"   ,rxsqr,
  23.    "TAN"    ,rxtan,
  24.    "TOPOWER",rxpower,
  25.    0        ,0
  26. };
  27.  
  28. /* Two utility functions not provided by the REXX calculator:
  29.    getdouble(argc)   checks that argc==1 and then gets that argument from the
  30.                      stack, returning it as a double.
  31.    stackdouble(argc) stacks the given double on the calculator stack and
  32.                      formats it according to the current settings */
  33.  
  34. double getdouble(argc)
  35. int argc;
  36. {
  37.    char *num;
  38.    char c;
  39.    int len;
  40.    double ans;
  41.    if(argc!=1)die(Ecall);
  42.    num=delete(&len);
  43.    if(len<0)die(Ecall);
  44.    
  45.    /* now remove all spaces (which are legal in REXX) for sscanf */
  46.    while(len&&num[0]==' ')num++,len--;
  47.    if(num[0]=='+'||num[0]=='-')
  48.       while(len>1&&num[1]==' ')num[1]=num[0],num++,len--;
  49.    while(len&&num[len-1]==' ')len--;
  50.  
  51.    /* Now the number is converted to double with sscanf */
  52.    num[len]=0;
  53.    if(sscanf(num,"%lf%c",&ans,&c)!=1)die(Enum);
  54.    if(ans==HUGE_VAL)die(Eoflow);
  55.    return ans;
  56. }
  57.  
  58. void stackdouble(num)
  59. double num;
  60. {  /* the number is simply sprintf-ed, stacked, and formatted. */
  61.    void rxformat();    /* the REXX format() function */
  62.    char result[25];
  63.    sprintf(result,"%.18G",num);
  64.    stack(result,strlen(result));
  65.    rxformat(1);
  66. }
  67.  
  68. /* The following functions implement the REXX math functions using
  69.    floating-point math. Some functions perform range checking first. */
  70.    
  71. int rxsin(name,argc)
  72. char *name;
  73. int argc;
  74. {
  75.    stackdouble(sin(getdouble(argc)));
  76.    return 1;
  77. }
  78. int rxcos(name,argc)
  79. char *name;
  80. int argc;
  81. {
  82.    stackdouble(cos(getdouble(argc)));
  83.    return 1;
  84. }
  85. int rxtan(name,argc)
  86. char *name;
  87. int argc;
  88. {
  89.    double arg=getdouble(argc);
  90.    double s=sin(arg);
  91.    double c=cos(arg);
  92.    if(c==0.)return -Edivide;
  93.    stackdouble(s/c);
  94.    return 1;
  95. }
  96. int rxasn(name,argc)
  97. char *name;
  98. int argc;
  99. {
  100.    double arg=getdouble(argc);
  101.    if(arg>1.||arg<-1.)return -Ecall;
  102.    stackdouble(asin(arg));
  103.    return 1;
  104. }
  105. int rxacs(name,argc)
  106. char *name;
  107. int argc;
  108. {
  109.    double arg=getdouble(argc);
  110.    if(arg>1.||arg<-1.)return -Ecall;
  111.    stackdouble(acos(arg));
  112.    return 1;
  113. }
  114. int rxatn(name,argc)
  115. char *name;
  116. int argc;
  117. {
  118.    stackdouble(atan(getdouble(argc)));
  119.    return 1;
  120. }
  121. int rxexp(name,argc)
  122. char *name;
  123. int argc;
  124. {
  125.    double val=exp(getdouble(argc));
  126.    if(val==HUGE_VAL)return -Eoflow;
  127.    stackdouble(val);
  128.    return 1;
  129. }
  130. int rxln(name,argc)
  131. char *name;
  132. int argc;
  133. {
  134.    double arg=getdouble(argc);
  135.    if(arg<=0.)return -Ecall;
  136.    arg=log(arg);
  137.    if(arg==HUGE_VAL)return -Eoflow;
  138.    stackdouble(arg);
  139.    return 1;
  140. }
  141. int rxpower(name,argc)
  142. char *name;
  143. int argc;
  144. {
  145.    double arg1,arg2;
  146.    if(argc!=2)return -Ecall;
  147.    arg2=getdouble(1);
  148.    arg1=getdouble(1);
  149.    if(arg2==0.){        /* Firstly, anything to the power 0 is 1. */
  150.       stack("1",1);
  151.       return 1;
  152.    }
  153.    if(arg1==0.&&arg2>0.){  /* Also, 0 to any positive power is 0. */
  154.       stack("0",1);
  155.       return 1;
  156.    }                    /* Otherwise return x**y = exp(y log x) */
  157.    if(arg1<0.)return -Ecall;  
  158.    arg1=exp(arg2*log(arg1));
  159.    if(arg1==HUGE_VAL)return -Eoflow;
  160.    stackdouble(arg1);
  161.    return 1;
  162. }
  163.  
  164. /* Two utility functions for the square root program below.  They provide
  165.    add and subtract functions which are quicker than going via REXX's
  166.    calculator. */
  167.    
  168. static int add(num1,len1,num2,len2)  /* add num 2 to num 1. */
  169. char *num1,*num2;                    /* Return the carry */
  170. int len1,len2;
  171. {
  172.    int i;
  173.    char c=0;
  174.    num1+=len1;
  175.    num2+=len2;
  176.    for(i=len2;i--;){
  177.       *--num1+=*--num2+c-'0';
  178.       if(c=(*num1>'9'))*num1-=10;
  179.    }
  180.    for(i=len1-len2;c&&i--;)
  181.       if(c=((++*--num1)>'9'))*num1-=10;
  182.    if(c)*--num1='1';
  183.    return c;
  184. }
  185.  
  186. static void sub(num1,len1,num2,len2)  /* subtract num2 from num1. */
  187. char *num1,*num2;
  188. int len1,len2;
  189. {
  190.    int i;
  191.    char c=0;
  192.    num1+=len1;
  193.    num2+=len2;
  194.    for(i=len2;i--;){
  195.       *--num1-=*--num2+c-'0';
  196.       if(c=(*num1<'0'))*num1+=10;
  197.    }
  198.    for(i=len1-len2;c&&i--;)
  199.      if(c=((--*--num1)<'0'))*num1+=10;
  200. }
  201.  
  202. int rxsqr2(name,argc)  /* Quick square root evaluator, for "low" precision */
  203. char *name;
  204. int argc;
  205. {
  206.    double arg=getdouble(argc);
  207.    if(arg<0.)return -Ecall;
  208.    stackdouble(sqrt(arg));
  209.    return 1;
  210. }
  211.  
  212. /* The following function implements the REXX square root function.
  213.    If the number of digits required is low enough, the floating-point
  214.    method (above) is used.  Otherwise, a digit-by-digit method is used. */
  215.    
  216. int rxsqr(name,argc)
  217. char *name;
  218. int argc;
  219. {
  220.    int n,m,e,l,z;
  221.    char x;
  222.    char *arg,*ans,*residue,*string,*tmp1,*tmp2;
  223.    int rlen,t1len,t2len;
  224.    extern int precision;  /* Extract the precision and the workspace from */
  225.    extern char *workptr;  /* the REXX calculator */
  226.    extern int eworkptr,worklen;
  227.    if(precision<=16)return rxsqr2(name,argc);
  228.    /* now the digit by digit algorithm for greater precision than FP math */
  229.    eworkptr=1;  /* Clear the work space, but allow for one digit of left-
  230.                    overflow (to be put at workptr[0]) */
  231.    if((n=num(&m,&e,&z,&l))<0)return -Enum;  
  232.    if(m)return -Ecall;
  233.    delete(&m);
  234.    if(z){       /* the square root of zero is zero. */
  235.       stack("0",1);
  236.       return 1;
  237.    }
  238.    /* workptr[n] is now the first in a sequence of digits of a positive
  239.       number.  The decimal point belongs after the first digit.  The
  240.       following code moves the decimal point to after the second digit,
  241.       and at the same time makes the exponent even. */
  242.    if(e%2)e--;
  243.    else workptr[--n]='0', /* (add leading 0) */
  244.         l++;
  245.    /* Expand the workspace to hold all the temporary material */
  246.    mtest(workptr,worklen,n+l+6*precision+20,n+l+6*precision+20-worklen);
  247.    arg=workptr+n;
  248.    if(l%2)arg[l++]='0'; /* make the length even by adding a trailing 0 */
  249.    ans=arg+l;
  250.    residue=ans+precision+1;
  251.    rlen=0;
  252.    string=residue+precision+precision+2;
  253.    tmp1=string+precision+2;
  254.    tmp2=tmp1+precision+3;
  255.    string[0]='0';
  256.    /* Now we have five numbers: ans, residue, string, tmp1, tmp2, all of
  257.       ample length, in the workspace as well as the argument. */
  258.       
  259.    for(n=0;n<=precision;n++){
  260.       /* invariants: length of arg used so far = 2n
  261.                      length of ans = n
  262.              length of string = n+1
  263.              length of residue = rlen <= 2n       */
  264.              
  265.       if(n+n<l)residue[rlen++]=arg[n+n],  /* Bring down two more digits */
  266.                residue[rlen++]=arg[n+n+1];/* from the argument to residue */
  267.       else residue[rlen]=residue[rlen+1]='0',
  268.            rlen+=2;
  269.        
  270.       /* find least x s.t. (string||x)*x>residue. Find (string||x-1)*(x-1). */
  271.       /* Invariant wil be: tmp1=(string||x-1)*(x-1) and tmp2=(string||x)*x  */
  272.       memcpy(tmp2,string,n+1);
  273.       tmp2[n+1]='1';
  274.       memset(tmp1,'0',t1len=t2len=n+2);
  275.       x='1';
  276.       /* Remove some leading zeros from the residue */
  277.       while(rlen>t2len&&residue[0]=='0')residue++,rlen--;
  278.       /* Below says "while(tmp2<=residue)" */
  279.       while(t2len<rlen||(t2len==rlen&&memcmp(tmp2,residue,t2len)<=0)){
  280.          memcpy(tmp1,tmp2,t1len=t2len);
  281.      string[n+1]=x;
  282.      if(add(tmp2,t2len,string,n+2))tmp2--,t2len++;
  283.      x++;
  284.      if(add(tmp2,t2len,&x,1))tmp2--,t2len++;
  285.       }
  286.       if(t2len>n+2)tmp2++; /* Restore tmp2 pointer to its original value */
  287.       
  288.       /* Append x-1 to the answer; subtract (string||x-1)*(x-1) from residue;
  289.          append 2(x-1) to the string. */
  290.       ans[n]=--x;
  291.       sub(residue,rlen,tmp1,t1len);
  292.       string[n+1]=x+x-'0';
  293.       if(x>='5')string[n+1]-=10,string[n]++;
  294.    }
  295.    stacknum(ans,n,e/2,0); /* We have the answer! */
  296.    return 1;
  297. }
  298.  
  299.